Ordered moment in the anisotropic and frustrated square lattice Heisenberg model 
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The two-dimensional frustrated next nearest neighbor Heisenberg model on the square lattice is a 
prime example for a spin system where quantum fluctuations can either destroy or stabilize magnetic 
order. The phase boundaries and staggered moment dependence on the frustration ratio J2/J1 of 
the exchange constants are fairly well understood both from approximate analytical and numerical 
methods. In this work we use exact diagonalization for finite clusters for an extensive investigation of 
the more general Jia,b-J2 model which includes a spatial exchange anisotropy between next-neighbor 
spins. We introduce a systematic way of tiling the square lattice and, for this low symmetry model, 
define a controlled procedure for the finite size scaling that is compatible with the possible magnetic 
phases. We obtain ground state energies, structure factors and ordered moments and compare with 
the results of spin wave calculations. We conclude that Jia,b exchange anisotropy strongly stabilizes 
the columnar antiferromagnetic phase for all frustration parameters, in particular in the region of 
the spin nematic phase of the isotropic model. 



PACS numbers; 75.10.Jm, 75.30.Cr, 75.30.Ds 
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I. INTRODUCTION 

The two-dimensional (2D) quantum J1-J2 Heisenberg 
model on a square lattice has recently found a number of 
realizations in layered V^^ {S = 1/2) compoundsir— . By 
using thermodynamic and high field^"^ investigations 
as well as NMR^'^, /xSR^, resonant X-ray scattering^ and 
neutron diffraction^i"— it was possible to locate these 
compounds in the phase diagram of this modeiii"— . It 
was concluded that all compounds are lying in the region 
of the columnar antiferromagnetic phase corresponding 
to ordering wave vector Q = (0, tt) or (tt, 0). 

The phase diagram is characterized by a single con- 
trol parameter, the frustration angle — tan^^(J2/«/i) 
and the energy scale is fixed by Jc = {Ji + J^)^ ■ The 
possible classical phases are a ferromagnet (FM) with 
TT — tan~^(l/2) < 4> < 37r/2, a Neel antiferromagnet 
(NAF) with -7r/2 < < tan-i(l/2) « O.IStt, and 
a columnar antiferromagnet (CAF) with tan^^(l/2) < 
(j> < TT — tan^^(l/2) « O.SSvr. In the transition re- 
gions NAF/CAF and CAF/FM frustration effects de- 
stroy the magnetic order in a small but finite inter- 
val. This can be concluded from exact diagonaliza- 
tion on small systemsi^Ti^l as well as spin wave calcu- 
lationsiSiiii^^iSl, series expansionS^Ti^S, and large- ex- 
pansion^ii^. It has been proposed that the true ground 
state in these regions is not a genuine spin liquid but 
exhibits hidden (nonmagnetic) order, namely a colum- 
nar dimerized phase with an excitation gap^^— and 
spin nematic order— respectively. The impact of spa- 
tial anisotropies on the columnar dimerized phase around 
(j> « O.IStt has been studied by exact diagonalization^l, 
the coupled-cluster method'^^ , and density- matrix renor- 
malization group methods^. The model also exhibits 
anomalies in frustration dependent magnetic and mag- 
netocaloric quantitie o"^^'^^ . 

The vanadates are in fact not strictly of square lattice 
type but small rectangular distortions lead to a small 



Jia — Jib anisotropy, however it was shown that it plays 
only a minor role for these compounds^. Furthermore the 
generalized anisotropic Jia,b-J2 Heisenberg model with 
large anisotropy has recently been invoked in the inter- 
pretation of spin wave excitations for the Fe-pnictide par- 
ent compounds—"—. It has also been discussed whether 
the observation of small ordered moments can be un- 
derstood within a frustrated local moment model. For 
a discussion of these topics we refer to Ref. H and the 
numerous references cited therein. However the discus- 
sion of ordered moment size in Ref. [l^ has so far been 
exclusively based on approximate analytical spin wave 
calculations. 

In this paper we want to endeavor a full scale numeri- 
cal approach based on exact diagonalization of finite size 
clusters to clarify the ordered moment size and its varia- 
tion with frustration and anisotropy effects in the Jia.b~J2 
Heisenberg model by an unbiased method. In addition to 
Jc and (f> this implies a further parameter 9 (which we will 
define later) characterizing the orthorhombic anisotropy. 

Our present work has two main objectives. Firstly 
it presents a new technical development how to apply 
the finite size scaling method systematically to a Heisen- 
berg model with low spatial symmetry. In most previous 
investigations of the isotropic J1-J2 model either restric- 
tive or not fully systematic methods have been chosen for 
the lattice tilings used for the scaling procedure^iii^ii^. 
In our case this has to be generalized because of the 
lower symmetry of the model and the appearance of non- 
degenerate phases with columnar magnetic order. We 
will discuss in detail how all possible tilings can be con- 
structed in a unique way, classify their symmetry and 
compatibility with classical ordered phases and introduce 
a precise concept of squareness to characterize their use- 
fulness for finite-size scaling analysis. In addition we will 
describe two different methods to calculate the staggered 
moment from the correlation functions and discuss their 
accuracy. 

Secondly, using these new technical developments we 
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will calculate the ground state energy, structure factor, 
and staggered moment as functions of frustration and 
anisotropy parameters and compare with classical behav- 
ior and the results from spin wave calculations. In par- 
ticular we will give a definite answer to which extent the 
ordered moment for general frustration and anisotropy 
parameters (0, 9) is modified as compared to the isotropic 
unfrustrated Neel state moment. We also discuss the in- 
fiuence of the anisotropy on the stability of the proposed 
spin nematic phase and show that a tiny deviation from 
the isotropic model reestablishes the columnar magnetic 
order. Further details on the physcial motivation for the 
investigations presented in this work are described exten- 
sively in Ref. [23- 

The paper is organized as follows: in Sect.|ll]thc model 
and its characteristic parameters are defined. In Sect. IIIll 
the technical implementation of the numerical analysis 
for the anisotropic model is described in detail. Then 
Sect. IIVI presents the finite-size scaling procedure, and 
Sect. |V] contains a discussion of the systematic depen- 
dence of the thermodynamic limit of the ground state 
energy Eq and the ordered moment Mq(Q) on the model 
parameters {4>,0). Finally Sect. IVll gives the summary 
and conclusions. 



II. THE MODEL 

The Hamiltonian for the two-dimensional = 1/2 
Jia,b-J2 model studied in this work is given by 

N N N 

(y>a 

The first two sums are taken over bonds between nearest- 
neighbor sites along the a and b directions of the rectan- 
gular lattice, respectively, and {{ij}) denotes bonds con- 
necting the next nearest neighbors along the diagonals of 
a rectangular plaquette. We use a parametrization of the 
exchange constants which facilitates the discussion of the 
whole phase diagram, and introduce a frustration angle 
(j) and an anisotropy angle 9. With these parameters, we 
define 

Jia — V2Jc cos cos 9, 
Jib — V^Jc cos (f) sin 9, (2) 
J2 = Jc sin (/i. 

Again Jc defines the overall energy scale of the model. 

The possible classical phases are ferromagnetic (FM, 
Qfm — (0,0)), Neel antiferromagnet (NAF, Qnaf — 
(tTjTt)) and columnar antiferromagnets (CAFa, QcAFa — 
(tt, 0), CAFb, QcAFb = (0,7r)). The latter are degen- 
erate for the isotropic {Jia — Ju) case with 9 = tt/A 
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FIG. 1. Classical phase diagram of the Hamiltonian defined 
in Eq. U]) as a function of the frustration angle (j> and the 
anisotropy parameter 6. For symmetry reasons, it is suffi- 
cient to restrict the discussion of the phase diagram to the 
parameter range — tt <<?!>< tt, 0<S<7r/4 indicated by the 
shaded area in the figure. 

or 9 ^ — 37r/4. Fig. [T] displays a sketch of the classi- 
cal phase diagram in the c^-^-plane. The shaded stripe 
in the plot denotes the parameter range — tt < ((> < tt, 
< 9 < tt/A which can be mapped onto the whole phase 
diagram applying discrete symmetry operations under 
which the Hamiltonian ([l} is invariant. 

Ref. f^i and references cited therein contain a thorough 
discussion of the classical phases and a spin- wave analysis 
of the ground-state properties of the Hamiltonian ((TJ . In 
particular, the inherent frustration leading to quantum 
fluctuations and possible moment reduction is investi- 
gated in detail. It is found that the staggered moment is 
stabilized in the columnar phases by introducing a spa- 
tial anisotropy, and a strong influence of frustration on 
the size of the moment can be excluded. However, linear 
spin-wave analysis a priori is a semiclassical theory, and 
naturally the question arises to what extent these results 
can be confirmed by strictly quantum-mechanical meth- 
ods. An unbiased but technically more involved approach 
is to determine the ground-state properties of the Hamil- 
tonian ([1]) by exact diagonalization on small finite tiles 
with periodic boundary conditions and to extrapolate the 
results to the thermodynamic limit using a scaling anal- 
ysis. 

III. TECHNICAL PREREQUISITES FOR A 
FINITE-SIZE SCALING ANALYSIS 

The most difficult aspect of finite size cluster calcu- 
lations is the implementation of an efficient finite size 
scaling procedure to obtain reliable values of physical 
quantities in the thermodynamic limit. In this section 
we present in detail the necessary ingredients. We de- 
scribe how we generate lattice tilings used for exact di- 
agonalization, classify them according to a newly intro- 
duced concept of their compactness or squareness (which 
will be defined later), their point- group symmetry, and 
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their compatibility with the classical phases of our model. 
These properties allow us to discriminate systematically 
the various clusters according to their usefulness for the 
finite size scaling analysis which is of central importance 
to obtain reliable results. Furthermore we will discuss 
the derivation of the set of wave vectors associated with 
a given tile and describe two different finite-size scaling 
procedures to calculate the ordered moment. 



A. Tiling the square lattice 

For a finite-size scaling analysis, it is common practice 
to select tiles according to certain geometrical or topo- 
logical properties. We first briefly describe three different 
schemes applied in the past to S = 1/2 spin models be- 
fore turning to the more general selection scheme used in 
this paper. 

Haan et alJ^ discuss the S — 1/2 nearest- neighbor 
Heisenberg antiferromagnet with helical boundary con- 
ditions and define an asymmetry parameter A — \£i — 
•^2 1/(^1 + -^2), where £i are the lengths of the edge vec- 
tors of the tiles under consideration. A square-shaped 
tile (considered as "good") has ^ = 0, but this is true 
for general diamond-shaped tiles, too. Those tiles having 
"small A" are selected for scaling. 

Restricting to strictly square-shaped tiles having at 
least C4 point-group symmetry is the recipe used by 
Schulz et alj^ for discussing ground-state energy and or- 
dered moment of the S* = 1/2 antiferromagnetic Ji — J2 
model. However, with this criterion only very few (two 
to four) tiles are eventually used for a linear least-squares 
fit. 

Another approach can be found in Ref. |45| . discussing 
the 5 = 1/2 nearest-neighbor XY and Heisenberg an- 
tiferromagnets with periodic boundary conditions. The 
authors introduce a parameter called the topological im- 
perfection of a tile, where a topologically perfect tile is 
defined as follows: A given lattice point on a tile con- 
tains ni nearest neighbors, n2 next-nearest neighbors, 
and so on, up to the imaxth-nearest neighbors where the 
sum over the rii reaches the tile area N. (Distance is 
measured as the minimal number of hops needed to get 
from one point to another.) If for all i < Jmax we have 
Hi = 4i, which is the number of ith- nearest neighbors on 
the infinite lattice, a tile is considered as topologically 
perfect. This concept is then extended to the notion 
of topologically perfect bipartite Neel lattices, i.e., the 
same conditions as described above are applied individ- 
ually to each of the two sublattices for antiferromagnetic 
Neel order. However, tiles are eventually chosen by hand 
in order to achieve a smooth finite-size scaling behavior 
of the ground-state energy per site and the square of the 
magnetization or staggered moment, respectively. 

The examples given above are in no way exhaustive, 
but illustrate one problem common to any finite-size scal- 
ing analysis, which becomes particularly apparent when 
applied to the full phase diagram of the Jia,b-J2 model. 



Firstly, only very few tilings might survive the final selec- 
tion, making a linear two-parameter fit to the ground- 
state energy or squared ordered moment of the Jia,b-J2 
model questionable, not to speak about higher-order cor- 
rection terms included in the fitting procedure.'*^ 

Secondly, the selection contains some arbitrariness 
which in our case would lead to selecting different tiles 
for scaling for different sets of exchange parameters, even 
within the same classical phase. 

Thirdly, given the edge vectors at of a particular par- 
allelogram, this is not a unique way to tile the square lat- 
tice. For example, upon replacing di by, say 2a2 — ai, we 
get a new tile with identical area which leads to an iden- 
tical structure of the resulting torus when introducing 
periodic boundary conditions. Even worse: In general, 
there are many possible generator matrices M = (ai, 02) 
of the same lattice tiling Am- 

However a more systematic approach to select tilings 
is possible. We will describe this scheme here and employ 
it for the Jia,b-J2 model. 

A basic requirement is a unique description of the lat- 
tice tilings To achieve this, we need some defini- 
tions. First, we introduce unimodular matrices U , which 
are defined as integer matrices of dimension s with de- 
terminant ±1. Let us further denote with [/('J) a unit 
matrix modified by having a single additional nonzero 
unit element at position (i,j), and introduce S*''^ as a 
unit matrix modified by having the (i,i)th element re- 
placed by —1. The unimodular matrices of dimension s 
form a (non-abelian) group Us under matrix multiplica- 
tion, and the generators of this group can be the elements 
and S'^'\ 

We need another definition; A s x s integer matrix 
H = (hij) is in Hermite Normal Form (HNF) if and only 
if 

hu>l, i = l .. .s, 
h^J =0, l<j<i<s, (3) 
h^j e [0, hu) , l<i<j<s. 

It can be shown*'' that (a) an arbitrary nonsingular inte- 
ger matrix M can be uniquely represented by a unimod- 
ular matrix U and an HNF matrix H such that 

U-M ^H, (4) 

and (b) that for two HNF matrices H and H' generating 
the lattice tilings Ah and Ah' , we have 

Ah = Ah'^H = H'. (5) 

In this way, we can uniquely represent an arbitrary lattice 
tiling of any given primitive Bravais lattice. 

The two-dimensional HNF matrices H have the form 

representing tiles with the special edge vectors hi = 
(/ill, hi2) and h2 = (0, /122) and area or number of sites 
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N = /iii/i22- For numerical purposes, it would already be 
sufficient to implement an algorithm using only Eqs. ([3]) 
and for building the list of possible tilings. 



B. Compact tiles 

In the spirit of Ref. [U where only square-shaped tiles 
have been used, we want to generalize this concept and 
utilize the notion of "squareness" or "compactness" of 
a tile for selecting the proper tiles for finite-size scaling. 
For this purpose, we introduce a parameter 




s/{s 



(7) 



for a nonsingular integer s x s matrix M, where Mi is 
the (non-square) matrix formed by dropping the ith row 
of M. p{M) measures the "compactness" ("squareness" 
in two spatial dimensions) of the s-dimensional paral- 
lelotope spanned by the row vectors of M: We have 
< p{M) < 1, and p(M) = 1 if and only if M de- 
scribes an s-dimensional cube, which we regard as the 
most compact lattice tiling in dimension s. 

However, calculating p{H) for the HNF representation 
of a lattice tiling is not very useful: According to Eq. (U), 
a single HNF matrix H represents a whole class Ch of 
tiles with matrix representation M which all lead to an 
identical lattice tiling Ku — Am- But in general, two 
matrices M ^ M' with Am = Am' have p{M) ^ p{M'). 
From each Ch, we therefore have to choose a tile which 
has the maximum compactness of all tiles of its class, 



Mc = M : p{M) 



max (p(Mr)) 



(8) 



and assign p{Mc) to the HNF representant H ~ UMc of 
its class Ch- 

With this scheme, we find 816 different classes of tiles 
with size N between 8 and 32. It is impossible to list all 
of them in this article. Instead, to illustrate the principle 
we display in Fig. [19] of Appendix [X] a list of all possible 
distinct ways to tile the square lattice using tiles of just 
the smallest useful size = 8. 

To label the tiles in a unique way, we use the scheme 
N:hii-hi2, where A^ is the number of sites or tile area 
and hij are the coordinates of the first edge vector of the 
HNF representation of a tile. The eight-site square for 
example has the label 8:2-2 (Fig. [191 third from below). 



C. Construction of the Brillouin zone for a finite 
lattice tiling 

To reduce the size of the matrices representing the 
Hamiltonian, it is useful to work with a basis reflect- 
ing the periodicity of the lattice, i. e., we construct Bloch 



states from the states of the {N, Sz) basis by applying 
Bloch's theorem. Algorithmically, the states of the origi- 
nal spinor product basis are collected into classes of wave 
functions which related to each other by translations with 
translation vector f and associated phase factor e'*^*". 

For the finite tilings we are using, the possible wave 
vectors k can assume only certain values. In this section, 
we discuss how to determine these wave vectors which are 
contained in the (first) Brillouin zone for a given lattice 
tiling. We use the fact that a translation by a reciprocal 
lattice vector will not change the phase of a wave func- 
tion, as required by Bloch's theorem. We construct the 
first Brillouin zone such that the origin is located in one 
of its corners. This is different from the commonly used 
Wigner-Seitz construction for infinite lattices, where the 
origin is in the center of the Brillouin zone. We cannot 
apply the Wigner-Seitz construction directly to an arbi- 
trary finite lattice, because we require the wave vector 
fc = to be part of the set of reciprocal lattice points 
generated, and the geometrical center of a given recip- 
rocal lattice tile constructed as described here does not 
necessarily have a wave vector associated with it. 

For the edge vectors oi and 02 of a tile and the corre- 
sponding basis vectors bi and 62 of the reciprocal lattice, 
we get from the orthogonality condition 

2n 



bi 













y an J 



(9) 



(10) 



022 
N \-a21, 
with integer coefficients aij, where 

A^= |Det (ai,a2)| 

is the number of sites or area of the tile under consid- 
eration. The reciprocal lattice vectors Gi — 2nex and 
G2 ~ 27re'y, where ex^y are the Cartesian unit vectors, 
are given by 

Gt = aubi + a2ib2, i = 1,2 (11) 

in terms of the reciprocal basis vectors. Within this co- 
ordinate system, the reciprocal lattice vectors Gi and G2 
span the parallelogram making up the first Brillouin zone, 
which contains exactly A^ wave vectors k — kibi + k2b2 
with integer coefficients ki- 

These wave vectors can be found with the following 
criterion: The projections of k onto the reciprocal lat- 
tice vectors Gi must be positive semidefinite (we put the 
origin (fci,A:2) = (0,0) into the lower left corner of the 
Brillouin zone) and less than the length of the latter, 
i.e., < k ■ Gi < Gi ■ Gi = A-k"^. With the expressions 
above for the basis vectors bi and the reciprocal lattice 
vectors Gi, we find as our defining relations for possible 
wave vectors k 

< Det(fc, G2) < A A < - Det(fc, Gi) < N, (12) 

writing the column vectors G; and k in the coordinate 
system spanned by the basis vectors bi- Conditions 
specify exactly A^ wave vectors k = fci6i -I- k2h2- Fig. [5| 
illustrates this procedure for tile number 14:1-11. 
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FIG. 2. Tile 14:1-11 in direct space (left), and the correspond- 
ing reciprocal lattice (right) according to Eqs. ([9| and (|lip . 
The two circles on the left-hand side mark the two equiva- 
lent points of the tile with maximum distance to the origin. 
Also shown is the first Brillouin zone of the lattice with edge 
vectors (reciprocal lattice vectors) Gi and G2. 




FIG. 3. Reciprocal lattice tile 14:1-11. The graphics shows 
the same tile as in the right-hand side of Fig. [2] corresponding 
to the first Brillouin zone, but here in the coordinate system 
spanned by the reciprocal basis vectors 61 and 62, indicated 
by the short arrows. The dotted circles denote the irreducible 
wedge of the Brillouin zone (tile 14:1-11 has C2 symmetry 
only), numbers give the size of the star of the corresponding 
wave vector. 



D. Space group symmetry 

To avoid unnecessary computations, we calculate 
Lanczos eigensystems only at those wave vectors k which 
are not related by a space group operation. Therefore we 
have to determine the irreducible wedge of the Brillouin 
zone. Figure 13] illustrates this for tile 14:1-11, which has 
C2 point group symmetry: The number of wave vectors 
is reduced from a total of 14 to eight nonequivalent ones. 
The small numbers in the figure denotes the size of the 
star of the corresponding wave vector. 

Since we are working with mappings of the full tile, it 
is useful to work in a coordinate system of the recipro- 
cal lattice vectors Gi and G2. For the projection of an 
arbitrary wave vector k — fciGi + k2G2 (not necessar- 
ily located in the first Brillouin zone) onto the reciprocal 
lattice vectors, we have 



k-Gi 1 

ki = = — Det 

IGiP N 



(fc,G2) 



(13) 



k. = 



k_G2 



-IdcI (k,Gi) 



(14) 



(compare Eq. (IT^ ). Note the minus sign in the last equa- 
tion and the interchange of Gi and G2: The projection of 
k onto Gi is proportional to the area of the parallelogram 
spanned by k and G2 and vice versa. 
We can write 



N 



Q<ni < N, mi eZ 



(15) 



for the coefficients of a wave vector k — M ■ q result- 
ing from the application of a space group operation M 
onto a wave vector (fin the first Brillouin zone. Mapping 
k back into the first Brillouin zone amounts to setting 
rrii = in the equation above. This gives us a convenient 
and numerically well defined procedure for mapping the 
Brillouin zone onto its irreducible part. 



E. Ordered ground states 

The classical Jia.b-J2 model on the square lattice has 
four ground states with ordering wave vectors Q = (0, 0), 
(tt, tt), and Q = (7r,0) or (0,7r). Although in the quan- 
tum case the corresponding wave functions, except for 
the ferromagnet, are not eigenstates of the Hamiltonian, 
it is important that the tilings of the infinite lattice are 
chosen such that these states corresponding to the clas- 
sically ordered phases are not suppressed when applying 
periodic boundary conditions. 

We can determine those tiles compatible with a clas- 
sical ground state by applying a test for the existence of 
the corresponding classical ordering vector Q in the list 
of wave vectors for a given tile. 



From Eqs. 
2n 



Tiibi + 71262 
, we get 



Hi e Z. 



ni 



0,22 
-0,21 



^^2 



/ ai2 
Van 



Qi 
Q2 



(16) 



(17) 



which has to be fulfilled for integer coefficients n^. 

• Ferromagnet: Qym — (0,0). Since we have chosen 
the phase of our wave functions such that fc = is 
always a valid wave vector, all tiles are compatible 
with the ferromagnet. 



CAFa phase: QcAFa = (tt, 0), and we have 

/ION 

ria = — . (18) 



ni 



Oil 

2 



Stated physically, the Mannheim distance (or 
Manhattan distance) c?m(pi,P2) = Sibii~-P2i|, 
counting the minimal number of hops on the lattice 
needed to get from one point to another, of any two 
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points Pi and pj belonging to the same sublattice 
projected onto the x direction of the lattice must 
be even, 



dxiVi^Pj) = \P3x -Pix\ = 2n, n€ Nq. 



(19) 



• CAFb phase: QcAFb — (0, tt). We get from condi- 
tion (HIl) 



ai2 0.22 



(20) 



so the components ai2 and 022 parallel to Cy must 
be even numbers. For completeness, the appropri- 
ate condition is given by 

dy{pi,pj) EE 2n, n E No (21) 

for any two lattice points pi and pj with 5f — Sj . 

• Neel phase: Qnaf = (tt, tt). Eq. ([TT]) can be solved 
with 

011+022 021-1-022 
"1 = 2 ' " 2 ■ ^ ' 

This is equivalent to 

dM{pi,Pj) = 2n, n e No 



(23) 



for any two lattice points Pi and pj on the same 
sublattice. 

• All four phases: All components of the edge vectors 
oi and 02 must be even individually in order to be 
compatible with the full set of classical phases. 

In this way, we can classify any given lattice tiling with 
respect to the classical phases it belongs to. 



For the C2V group, two sets of mirror "planes" exists: 
The point group C4V contains two isomorphic subgroups 
C2yroct and C2V'dia, corresponding to a tile with either 
a rectangular shape (mirror planes parallel to the edges) 
or a diamond-like shape (mirror planes along the diago- 
nals). This distinction is important when discussing or- 
thorhombic or trigonal symmetry breaking, since C4Y is 
reduced to the respective C2V symmetry in this case. For 
simplicity, we also label those tiles having C2V but not 
C4V symmetry accordingly in Table HI 



G. Static structure factor and ordered moment 

In a finite system at zero magnetic field the ground 
state has Sz = and therefore does not exhibit the spon- 
taneous symmetry breaking of the infinite lattice. The 
ordered moment M{Q) rather has to be obtained indi- 
rectly from the properly normalized static structure fac- 
tor according to 



SNiQ) 



1 ^ 



N 



(24) 



N 



i=2 



where the angular brackets denote the ground-state ex- 
pectation value in the Sz = 0, k — subspace of the 
Hilbert space. In the thermodynamic limit, if Q is the 
ordering vector of the corresponding classical phase, we 
can then identify 

M^iQ)^ lini Af|r(Q) = C(Q) 1™ Sn{Q) (25) 



Selection of tiles 



with the appropriate normalization Af of Sn{Q)- Here 
we have introduced a factor 



In order to discuss spatial anisotropics, we do not re- 
gard lattice tilings as equivalent which are related by a 
point-group operation on the square lattice. An exam- 
ple would be the four different tiles 8:1-2, 8:1-6, 8:2-1, 
and 8:2-3, see Fig. [T21 In this way, we get 816 differ- 
ent tilings of the square lattice with tiles between 8 and 
32 sites. Out of these, we select, for each even tile area 
and for each classical phase, the tile having the maxi- 
mum squareness. In addition, for N = A£, i € N, we in- 
clude the tile with maximum squareness containing both 
CAFa and CAFb ordering vectors, which is equivalent 
to containing all four classical ordering vectors. These 
special tiles are of particular importance for the spatially 
isotropic model, due to the degeneracy of the two colum- 
nar phases in this case. The resulting list is displayed 
in Table HI Each line contains the tile label, its compat- 
ibility with classical phases (Sect. IIII E[) . its squareness 
fSect. HITB]) . and its point group symmetry (Sect. UlTDt . 



C(Q) 



1, Q = or (tt, tt) 

2, Q = (^,0) or (0,7r) 



(26) 



to account for the additional lattice rotation symmetry 
breaking in the CAFa and CAFb phases. 

For a perfectly ordered classical state with wave vec- 
tor Q, the ordered moment assumes its maximum value, 
M{Q) = S, and we have 

(5,5^,) = =52g-iQ(fl.-i?,)^ (27) 

Assuming perfect order for the finite tile under consider- 
ation, too, we have 



Sn{Q) 



N 

Jf 
1 

Jf 



[S{S+1) + {N -1)S^ 
NS {NS + 1) . 



(28) 
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Tile NAF CAFa CAFb □ C2 C2Vroct C2Vdia d dv 



oil- J, 


• • 




1 nrtn ~ 
1.000 • 


• 


• 


• 


• 


10:1-3 


• — 


— 


1.000 • 






• 




10:1-4 




* 


0.966 • 










1 n.o 
iU:z-z 






n nan ~ 
0.966 • 










12:3-0 




* 


0.980 • 


• 








12:4-0 






0.980 • 


• 








12:1-5 


• — 




0.960 • 




• 






10.0 


• • 




n oni . 










1/1.1 


• — 




0.961 • 










1/1.1 /I 
i4:i-4 




* 


0.938 • 










1 /1 .0 Q 






U.yoo • 










lo:4i:-U 






i.UUU • 


• 


• 


• 


• 


18:3-3 


• — 




1.000 • 


• 


• 


• 


• 


10.1 /I 




* 


0.9/0 • 










ICO /I 






u.y/o • 










or».o /i 


• • 




1.000 • 






• 




22:1-6 




— 

* 


0.981 • 










22:2-4 






0.981 • 










22:1-5 


• — 




0.961 • 










24:1-10 




* 


0.988 • 










24:2-5 






0.988 • 










/I . 1 T 


• — 




n non — 
0.980 • 




• 






/I . /I r» 


• • 




A nc;n — 
0.960 • 


• 








26:1-5 


• — 


— 


1.000 • 






• 




o^^. 1 in 






n nfi /I . 










zd:2-o 






n /I . 
U.964 • 










28:1-8 




• 


0.986 • 










28:4-3 






0.986 • 










28:2-4 


• • 




0.961 • 










30:5-0 






0.992 • 


• 








30:6-0 






0.992 • 


• 








30:1-5 


• — 




0.974 • 










32:4-4 


• • 


• 


1.000 • 


• 


• 


• 


• 



TABLE I. Label, classical phase compatibility, squareness, and point groups for selected lattice tilings between 8 and 32 sites. 
For each even area A'' and for each classical phase, the list contains the compatible tile with maximum squareness as defined in 
Eq. ((7|). For A'' = 12 and 24, the tiles compatible with all classical phases are included, too, although they have a comparatively 
small squareness. The tile labels have the form N:hii-hi2, where hij are the components of the first edge vector of a tile in 
HNF representation. Those tiles compatible with all four classical phases, required for the discussion of the spatially isotropic 
model with columnar order, are underlined and typeset in bold. 



If we require Sn{Q) — also in this case, we have to 
set 



(29) 



which is the normalization we use for any tile included 
into our finite-size scaling_analysis for Mn{Q). This is 
in accordance with Refs. 21:1 andlHll and slightly deviates 
from the J\f = N'^ normalization commonly used by many 
authors. 



In the FM regime, although the fully polarized (all-up) 
state is an eigenstate of the Hamiltonian, the structure 
factor at the antiferromagnetic ordering vectors remains 
small, but finite. Assuming perfect order again, we get 
for the individual terms in the sum Eq. p4p 



SiSi) ^S{S + 1), (SiS,) = S', j^l, (30) 
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leading to 



Sn{Q) 



N 



N 

5(5+l) + 52^e- 
i=2 



(31) 



for all three antiferromagnetic wave vectors in the ferro- 
magnetic phase. As required, this value vanishes in the 
thermodynamic limit. 



H. Long-distance correlations 



Let us restrict to tiles with an even area N which is a 
necessary condition for being compatible with at least 
one of the non-FM phases of the Jia,b-J2 model. The sum 
in the above equation contains only exponentials which 
can acquire the values +1 or —1 for Q = (tt, tt), (tt, 0), or 
(0, tt). For each of these three wave vectors, there are N/2 
sites j with distance Rj to site 1 which have phase +1, 
and N/2 sites with phase —1. Site 1 obviously belongs to 
the former, such that we have N/2 — 1 terms left in the 
sum over the exponentials above having phase +1, and 
we get 



^ = ( ^ _ 1 j X (+1) + ^ X (-1) = -1. (32) 



With J\f given by Eq. (P^)) , we therefore have 



SNiQ) 



S 



TV + 1/5 



(33) 



An alternative way of calculating the ordered moment 
in the thermodynamic limit is mentioned for example 
m Refs. ^ and ^sl For the infinite system, in an or- 
dered phase with a staggered moment, the spin correla- 



tion function factorizes for large distances \Ri - 
Eq. ([27]) simplifies: 



Ri\, and 



lim 



S-i S 'j 



5, 



Af2(Q). (34) 



Working with a finite compact tile, we can extrapolate 
the spin correlation function for a single pair of spins, 
defining lattice points i and j such that their distance 
on the tile is maximized under given periodic boundary 
conditions. Without loss of generality, we can restrict 
ourselves to finding the pair (l,j) or just site j being 
maximally apart from the origin. 

We have to define a metric on a tile refiecting the pe- 
riodic boundary conditions: Each tile has four corners 
with coordinates (0, 0), ai, 02, and ai -I- 02, which are all 
equivalent to point 1. Using the Mannheim distance cZm 
again, we define the toroidal Mannheim distance between 
two points Ri and Rj on a tile as 



-{i,j) = minjdM (^Rij,Oj , (-Rij,ai) , c^m (^Rij,d2^ ,du 



Rii, fll + a2 



)} 



(35) 



yielding the minimal distance between two points Rt and 
Rj on a torus. Here the vector Rij is the point emerging 
from Rj — Ri when shifting Ri to the origin, projected 
back into the original tile: With 

Rj - Ri ^ {xji - xa)ai + {xj2 - Xt2) 0,2 

= ]^ [("ji " ai + ("-J2 - ni2) 02] (36) 

this point is simply given by 

Ri] = [((%i - f^il) mod TV) ai (37) 
+ {{nj2 - ni2) mod Af) 02] . 
With respect to the distance dx, we can then define 



j"=max<j^(i_3)0"gT(Af)) 



(38) 



For a square with size N = L'^, L even, the point Rj = 
(L/2,L/2), i.e., the center of the square, has maximum 



distance dx from the origin. However, in most cases, due 
to the lack of a lattice point located in the geometrical 
center of the tile, there will be, for even N, two sites j 
having the same maximum distance dx from the origin, 
see Fig. [2] for an examplei^ For our calculations, we just 
select one of them. We then can give an estimate for the 
ordered moment as 



M^(Q) = lim 



mUQ)- 



(39) 



IV. FINITE-SIZE SCALING ANALYSIS 

In Refs. [m and [H, the area dependence of ground- 
state properties of the two-dimensional antiferromagnetic 
Heisenberg model has been derived using chiral pertur- 
bation theory. In particular, for the ground-state energy 
and the ordered moment, the following scaling behavior 
has been found: 



EqN — Eq 



1 



7V3/2 



1 

4piV2 



1 



,(40) 
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Ml{Q)^M^{Q) + a 



cx± iVi/2 



O 



(41) 



where c = \fpJxL is the spin-wave velocity, p the spin 
stiffness constant, and XJ- the transverse susceptibihty. 
a = 0.620704 and (i = -1.437745 are numerical con- 
stants. 

Although these results were proposed only for the 
isotropic unfrustrated case, we use them for the whole 
phase diagram. We argue that the form of the scaling 
functions should not depend on range and anisotropy of 
interactions within certain limits, as long as the model 
belongs to the same universality class. However the in- 
dividual coefficients in Eqs. (|40|) and (|41|) might change. 
Hence we assume the following size dependences for the 
ground-state energy and the ordered moment: 



'ON 



Eo 



62 



Ml{Q) = AP{Q) 



iVl/2 



m| 
N ' 



(42) 
(43) 



where the latter scaling function is applied to both 
Mm{Q) and Mn{Q)- With the exception of the point 
{4>,9) = (0,7r/4) (isotropic nearest- neighbor exchange Ji 
only), for all combinations of </> and 6 discussed here we 
use only the first two terms in the equations above for our 
scaling analysis. This is due to the comparatively small 
number of different areas of the tiles actually included in 
the calculation. 

We have calculated the ground-state energies Eon^ 
structure factors Mf^{Q), and long-distance correlation 
functions M^{Q) at the respective ordering vectors for 
tilings of 11 different sizes N between 8 and 28 at more 
than 300 points in the (0, 9) plane, producing roughly 
17 500 data sets altogether. To demonstrate the proce- 
dure, in the following we show a few examples at some 
significant points in the phase diagram before presenting 
our main results in section IVl 



A. The ground-state energy 

The ground-state energy of the Hamiltonian ^ is cal- 
culated in the subspace with total spin 5^ = at wave 
vector k = Q (the respective classical ordering vector) 
except for tiles with area N = 4£, £ G N, where the 
ground state is found in the k = sector. Figs. HHH show 
the dependence of the ground-state energy per site Eon 
on the reciprocal area of the tiles for selected model pa- 
rameters {(j), 9). (We plot —Eqn in order to be consistent 
with previous studies^i^.) The insets display schemati- 
cally the position in the phase diagram the plots refer to. 
Here and in the following, energies are measured in units 
of Jc, unless mentioned otherwise. 

The short horizontal dashes show the calculated en- 
ergies for those tiles compatible with the corresponding 
classical phases for a given parameter set {(j), 9). The en- 
ergy eigenvalues of those tiles having maximum square- 



ness are indicated by the small open circles. We have 
also determined Eqn for tiles incompatible with the clas- 
sical phases. In these cases, we get consistently much 
higher values for Eon, and we omit them in the plots. 
The larger value of the energy for these tilings is due to 
the fact that the interactions governed by the geometry 
of the tile prevent achieving a lower value of the energy. 

Furthermore there are tiles which are compatible with 
classical phases but correspond to very skewed parallel- 
ograms with small squareness (less than 1/2), in some 
cases they are ladders or even chains. They have a much 
lower ground-state energy (shown as dashes above the 
scaling line in Figs. |4] and [5]) than the tiles of equal size 
and larger squareness. Although the energies of these 
finite clusters with low squareness are shown in the fig- 
ure for completeness they are not used for the scaling 
procedure. 

For each tile area iV, we have to choose one value for 
Eon out of the stack of eigenvalues to be included into 
the finite-size scaling fit applying Eq. p2)) . According 
to our selection criteria, we choose the tile which (i) is 
compatible with the classical phase the parameter set 
(0, 9) belongs to, and (ii) is the best approximation to a 
square in the sense of Eq. ([7]). See Table U for a list of 
them. The values Eqn obtained in this way, labeled by 
the open circles in Figs. HHH are then used for a fit to 
Eq. P^ . indicated by the solid lines in the figures. 



1, Isotropic nearest-neighbor exchange 

The top of Fig. |4] shows the ground-state energy scal- 
ing of the conventional unfrustrated antiferromagnetic 
Heisenberg model with Jia = Ju = Jc and J2 = 0, 
equivalent to {(p,9) — (0,7r/4). We obtain the value 
Eo — —0.66(5) for the extrapolated ground state en- 
ergy which is in agreement with previous calculations, 
see, e.g., Refs. Hi] and lisi and references cited therein. 
It should be noted that not always the tile having the 
highest squareness also has the highest ground-state en- 
ergy compared to the other tiles with the same area N. 
However, the differences are small, and not visible on the 
scale of Fig. |4l 



2. Next-nearest neighbors and spatial anisotropy 

The lower plot in Fig.|3]shows the scaling of the ground 
state energy for the point (0,0) = (— 0.37r, 7r/4) corre- 
sponding to the isotropic model with finite ferromagnetic 
next-nearest-neighbor exchange J2. Ferromagnetic J2 
stabilizes the (tt, tt) order, and we obtain a lower ground 
state energy Eo — —0.76(2) accordingly. 

The top of Fig. [5] illustrates the scaling behavior in 
the CAF phase of the isotropic model with {(j),9) = 
(0.37r, 7r/4). Here all exchange constants are antiferro- 
magnetic, and the competition between them leads to 
columnar order. Because Jia = Ju here, we have two 
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FIG. 4. Ground-state energy of the isotropic model as a func- 
tion of Top: Conventional Neel AF Heisenberg model 
{Jia = Jib, J2 = 0). Bottom: Additional ferromagnetic next- 
nearest neighbor exchange J2 < 0. The dashes indicate the 
energy per site for different tilings. Those having maximum 
squareness, marked with a circle, are used for the scaling 
and extrapolation. The solid line shows the scaling curve, 
Eq. (|42p . The insets in these and the following plots show 
the parameter position {<j), 9) in the phase diagram, cf. Fig. [1] 
The numbers near the abscissae denote the tile size A'^. 



equivalent ordering vectors (tt, 0) and (0,7r), labeled by 
CAFa and CAFb, respectively. Only the tiles which con- 
tain both the CAFa and CAFb classical phases have truly 
the symmetries required by the Hamiltonian in this case. 
Consequently only the tiles with size N = Al, I € 'H 
are acceptable for scaling, and we have much less data 
points available as compared to the Neel phase. Fig. S) 
As before, we select out of these the tiles with maxi- 
mum squareness for the scaling procedure which leads to 
a ground state energy Eq = —0.53(2). 

In the anisotropic case, CAFa and CAFb phases are 
no longer equivalent. An example of CAFa for the max- 
imally anisotropic case (</>, ^) = (0.27r,0), where Ju = 
and J2 antiferromagnetic, is shown in the bottom part 
of Fig. m Here the tiles do not have to be compatible 
with the (0, tt) ordering of the CAFb phase, therefore 



FIG. 5. Finite-size scaling of the ground-state energy in (top) 
the CAF phase of the isotropic and (bottom) the CAFa phase 
of the maximally anisotropic model for different tilings. La- 
bels, symbols and lines have the same meaning as in Fig. [l] 



the number of possible tilings is larger. Also the scaling 
behavior is improved again. Compared to the isotropic 
case, the extrapolated value for the ground-state energy 
is lower, indicating that the introduction of a rectangular 
anisotropy stabilizes the columnar order. 



3. Magnetically disordered regimes 

The scenarios discussed until here have one property 
in common: All of them have parameter sets (0, 9) which 
are located deeply inside the corresponding classically 
ordered phases, and all of them show a good and well- 
defined scaling behavior of the ground-state energy. (We 
give a precise definition and discussion of the meaning of 
"good" and "well-defined" in Sect. HVCl ) 

However, this changes when approaching the disor- 
dered regimes at the borders of the columnar phases. 
Even with our stringent conditions for tile selection, the 
behavior in the region where the transition from the Neel 
to the CAF phase occurs is quite different. In fact an area 
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FIG. 6. The tile size dependence of the ground-state energy 
in the disordered regime at the CAF-FM border. In this case 
a reUable Unear scaUng result cannot be obtained from the 
circles. 



dependence as given by Eq. no longer seems to apply, 
and the concept of choosing the most-square-like tiles for 
scaling apparently becomes inappropriate. An example 
of this behavior is displayed in Fig. [6] corresponding to 
the well know disordered case of the isotropic model with 
J2/J1 « 1/2. 

Expressed differently, finite-size scaling breaks down 
and is no longer a useful concept by itself in and near the 
disordered regions at the edges of the columnar phases 
in the phase diagram. This will be described in more 
detail including a discussion of the correlation functions 
and the ordered moment in section fVl 

We conclude that a stable finite-size scaling analysis of 
the ground state energies can be achieved inside the NAF 
and CAFa,b regions for all frustration and anisotropy pa- 
rameters provided a careful selection of tiles according 
to phase compatibility and maximum squareness in each 
case is carried through. However the scaling analysis can- 
not be made in the transition regions close to the phase 
boundaries in Fig. [T] This is correlated with the break- 
down of the ordered moment in these regions as will be 
shown in the next section. 
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FIG. 7. Finite-size scaling of the ordered moment for the 
isotropic nearest-neighbor AF Heisenberg model for differ- 
ent tilings. Top: M|r (Qnaf) derived from the structure fac- 
tor 5'jv(Qnaf)- Bottom: M^(Qnaf) derived from the long- 
distance correlation function. The horizontal dashes and open 
circles denote the values for different tiles. The latter are 
those with maximum squareness which are used for the scal- 
ing and extrapolation. The solid lines show the fit of Eq. (|43p 
to M^(Qnaf), the dashed lines the fit to M'^{QjqAF)- To fa- 
cilitate the comparison, we reproduce both fits in both plots. 



Again, apart from the conventional Neel antiferromagnet 
we apply a linear scaling only, ignoring the last term in 
Eq. dMl). 



B. Ordered moment, structure factor and 
correlation function 

Using the ground-state wave-function obtained by di- 
agonalizing the Hamiltonian matrix, we calculate the 
expectation values of the spin correlation functions 
(SiSj). Subsequently, we calculate the finite-size equiv- 
alent Mff{Q) of the ordered moment applying the two 
methods described in Sects. HTl Gl and UlI HI In the same 
way as for the ground-state energy in the previous sec- 
tion, we select the tile with maximum squareness for each 
tile area N and fit Eq. (^5)) to the resulting data points. 



1. Isotropic nearest-neighbor exchange 

The top part of Fig. [7] shows the ordered moment de- 
rived from the structure factor, Sect. IIII G[ the bottom 
part shows (Q) obtained from the long-distance corre- 
lation function, Sect. IIII Hi for the unfrustrated (J2 = 0) 
isotropic AF Heisenberg model {(f> — 0,6 = 7r/4). As be- 
fore, the horizontal dashes show the values for different 
tilings, and the tiles with maximum squareness are indi- 
cated by circles. From the top figure it is evident that 
the most square-like tiles have also the largest structure 
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factor at the ordering vector Q — (tt, tt) for the NAF 
phase. 

A fit of Eq. (|43p using all three terms has been ap- 
plied to both M^{Q) and M^{Q) separately. The solid 
lines in the two plots of Fig. [7] denote the result of fitting 
to Mff{Q), and we get M(Qnaf) = 0.30(3). This is in 
good agreement with previous studies^iii^i^. The same 
fit applied to M'^{Q), indicated by the dashed lines in the 

figure, gives a slightly lower value M(Qnaf) = 0.27(6) 
for the thermodynamic limit. We believe this to be a 
good indicator for the accuracy within which we can de- 
termine asymptotic values for the ordered moment. Also 
the relative error (to be explained below) is larger for 
the latter method, which is a consequence of the fact 
that only a single correlation function is evaluated, while 
in the first method, a Fourier transform using all possible 
(SiSj) is taken. 



ors and spatial anisotropy 



2. Next-nearest 



The ordered moment scaling for ferromagnetic J2 < 
{4> — — 0.37r), both in the isotropic (Jia = Jib) and max- 
imally anisotropic {Jib = 0) case is shown in the top 
and bottom of Fig. [S] respectively. The horizontal dashes 
denote M'^{Q) for individual tilings, the solid line rep- 
resents a fit with Eq. (|43p to the values for most-square- 
like-tiles (circles in the plot). The dashed lines denote fits 
to M'^{Q) for the same set of tiles. (Individual values are 
not shown.) Again, a comparison of the extrapolated val- 
ues for M^{Q) from the two different scaling procedures 
can serve as an indicator of the quality of the finite-size 
scaling analysis. 

Fig. |9] shows the scaling of the ordered moment in 
the columnar phase for antiferromagnetic J2. For the 
isotropic case (top), as before only those tiles compatible 
with both columnar phases, equivalent to compatibility 
with all four classical phases can be used. In the infinite 
system, which has C^v point-group symmetry, the wave 
vectors (tt, 0) and (0,7r) are equivalent. However, most 
finite tiles have a spatial symmetry lower than C4 (see 
Table [l|, meaning that the equivalence between (tt, 0) 
and (0, tt) is lost. For this reason we take the sum of the 
structure factor at these two wave vectors and use the re- 
sulting data points for scaling. This is not necessary for 
the anisotropic case {9 ^ 7r/4, — 37r/4), of which an ex- 
ample is shown at the bottom of Fig. [HI since then CAFa 
and CAFb are different phases anyway. Here also a larger 
number of clusters is available for the scaling. Note that 
in Fig. ini we plot the area dependence of the structure 
factor and the long-distance correlation function includ- 
ing the factor CiQ) = 2 introduced in Eq. (pS)) . which 
is due to the spatial symmetry breaking induced by the 
columnar order. 
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FIG. 8. Finite-size scaling of the ordered moment in the NAF 
phase with ferromagnetic J2. Top: isotropic model, Jia = 
Jib- Bottom: maximally anisotropic case, Jit = 0. The 
horizontal dashes denote M%{Q) for different tilings, the solid 
lines are fits with Eq. (|43p to the values for most-square-like- 
tiles, marked with a circle. Results for M'^{Q) for individual 
tiles are not shown, fits to Mff{Q) are displayed as the dashed 
lines. 



3. Magnetically disordered regimes 

As it has been the case for the ground-state energy, in 
the disordered region between the Neel and the colum- 
nar phase of the isotropic model (Jia = Jib — Ji, 
J2/J1 ~ 1/2), the behavior of both the structure factor 
and the correlation function do not show a systematic 
dependence on the tile size. This is obvious from Fig. [TUl 
which shows the calculated values for both M^{Q) (top 
part) and Mff{Q) (bottom part). Disregarding for once 
the squareness of our tiles, it appears that there are two 
"stripes" of values as a function of 1/Vn which both 
seem to extrapolate roughly to 0. However, a quantita- 
tive analysis in the sense of Eq. cannot be made. 

We observe a very similar behavior also for J2/J1 ~ 
— 1/2 (with ferromagnetic Ji) in the spin-nematic re- 
gion of the phase diagram. In order to see more details, 
the long-distance correlation function (SiSm), used in 
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FIG. 9. Finite-size scaling of the ordered moment in the 
columnar phase with antiferromagnetic J2. Top: isotropic 
model, Jia = Jib- Bottom: maximally anisotropic case, 
Jib = 0. The horizontal dashes denote M%{Q) for differ- 
ent tilings, the sohd lines are fits with Eq. (|43l) to the values 
for most-square-like-tiles, marked with a circle. Results for 
M^{Q) for individual tiles are not shown, fits to M%{Q) are 
displayed as the dashed lines. The extrapolated values ob- 
tained from the scaling correspond to half the value of the 
ordered moment, due to the spatial symmetry breaking intro- 
duced with columnar ordering. 



Eq. ([55]) . in this narrow region between CAF and FM 
phase is plotted as a function of the frustration angle 
in Fig. [11] for tiles of different size. The labels of the 
individual tiles are given in the legend of the plot. The 
solid vertical line indicates the classical phase boundary. 
Apart from the eight-site tile, which is too small anyway 
because sites 1 and m are just two hops apart, for each 
tile (SiSm) is monotonously decreasing before jumping 
to the ferromagnetic value (SiSm) ^ = 1/4. Those 
correlation functions which are ferromagnetic (positive, 
both sites 1 and m on the same sublattice) in the colum- 
nar phase even change their sign before the jump. This 
sign change possibly can serve as an indicator of the 
breakdown of columnar order. However, the transition to 
ferromagnetic correlations in the region of 0.881 tt < < 
0.89l7r seen here is not happening uniformly for all tiles. 



FIG. 10. Tile area dependence of the ordered moment in 
the disordered regime at the NAF-CAF border for isotropic 
nearest-neighbor exchange. The values are obtained from 
(top) the structure factor and (bottom) the long-distance cor- 
relation function, and the circles show the values for the tiles 
with maximum squareness. No finite-size scaling is possible. 
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FIG. 11. Dependence of the long-distance correlation function 
on the frustration angle <j) in the spin nematic region between 
the CAF and FM phases of the isotropic model for those tiles 
with maximum squareness which are compatible to all four 
classical phases. The legend in the plot lists the individual 
tiles and their symbols. 
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hence it cannot be used to extrapolate to the thermody- 
namic hmit. This wih be discussed more quantitatively 
in the next section. 

We can draw the same conclusions for the ordered mo- 
ment scaling as we did before for the ground state energy: 
If we carefully select the allowed clusters with appropri- 
ate classical phase compatibility and maximum square- 
ness we can carry out a well defined scaling procedure to 
the thermodynamic limit for the stable NAF and CAFa,b 
regions. However in the disordered regions at the borders 
of the CAFa/b honeycomb (see Fig. [Ij scaling is impos- 
sible, which is an indication of the imminent breakdown 
of magnetic order due to quantum fluctuations. Similar 
conclusions were obtained already from linear spin wave 
theory^!. 



C. Accuracy of extrapolated values 



ein = 0.25 



For the linear scaling (Eqs. (j42|) and (|43l) without the 
last term) applied to extrapolating the results to the in- 
finite lattice, an analysis of the quality of the scaling is 
essential. There are several measures that can be in- 
troduced to characterize the quality of the scaling. A 
convenient measure is the relative error in the (xi-norm, 
generally defined as 



ll/lloo 



10" 



(44) 



||{a;i,a:2, • • -llloo = max(|a;i|, \x2\, ■ ■ ■) ■ 



Here, / is the set of fitted values at the points 1/A^, and 
d the set of data calculated for the maximum-squareness 
tiles with area N . The exponent p in the above equation 
can be interpreted as the number of significant digits^ 
for the extrapolated value fo at iV — !• oo. 

Figure [H] shows a plot of Crei for Mfj{Q) of the 
isotropic model as function of the frustration angle 0. In 
the ferromagnetic phase, the fully polarized state is the 
ground-state and the error reflects the numerical noise, 
i. e., the accuracy of the floating point operations, which 
are of order 10~^^, and has been excluded from the fig- 
ure. (We have determined M^(Q) for (0, 6) values in the 
ferromagnetic region mainly in order to verify the correct- 
ness of our numerical implementation of Eq. (1431) .) 

In the well-ordered region of NAF and CAF we have 
Eroi — C'(10^'^...10^^), and we can regard the scaling 
procedure in these areas as stable. This is in contrast to 
the magnetically disordered regions at the borders of the 
columnar phases, where the strong increase in erei clearly 
indicates that the scatter of the points is much higher 
and makes the extrapolated value unreliable. The solid 
horizonal line in Fig. [T^ denotes ej-oi = 0.1. According to 
Eq. (|44p . this is the maximum error at which the extrapo- 
lated values for ]VP{Q) have at least one significant digit. 
In other words, a magnetic order parameter M'^{Q) to be 
used to characterize the nature of the ground state can- 
not be obtained anymore in those regions with e^ci > 0.1. 




0.001 



FIG. 12. Relative error e^ei as defined in Eq. (1441) of the fits to 
Mjf (Q) for the isotropic model. Note the logarithmic scale of 
the ordinate. The error in the magnetically disordered regions 
is at least one order of magnitude larger than in the ordered 
sectors of the phase diagram. The solid horizontal line denotes 
the value erei = 0.1, indicating the maximum error acceptable 
for having at least one significant digit in M^{Q). 



-0.2 
-0.3 
-0.4 
^-0.5 
-0.6 
-0.7 
-0.8 



;/;!■ = 0.25 



FM 


NAF 


CAF 


FM: 


- Classical 
















— LSW 








-i 








ED 



























-1.0 



-0.5 



0.0 
cf>/n 



0.5 



1.0 



FIG. 13. The ground-state energy as function of the frustra- 
tion angle (p for the isotropic model with fixed 6 = n/A . The 
classical energy is shown as dashed line, and the spin-wave re- 
sults including zero-point fluctuations are presented as solid 
line. Dots indicate the values for the energy Eq obtained from 
extrapolating our exact-diagonalization data. 



V. RESULTS 

Here we describe the systematic dependence of the 
thermodynamic limit of the ground state energy Eg and 
the ordered moment Mq{Q) on the model parameters 
((/), 9) as obtained from our finite-size analysis discussed 
in the previous section. 



A. Ground-state energy 

Fig. [T3] shows the energy dependence on the frustra- 
tion angle for the isotropic model with O/tt = 1/4. 



15 



(j>l7l = 



-0.2 




FM 


CAFa 


NAF 


CAFb 


-0.3 






-< 


y- 


-< 








-0.4 




y- 


-< 


y- 


4 / 


\ 


\ 


















\ 












-0.5 








V'- 






-0.6 




- Classical 










T 










-0.7 




ED 








-0.8 















-1.0 -0.5 0.0 0.5 1.0 

e/n 



FIG. 14. Ground state energy for the anisotropic model with 
first neighbor interaction only (J2 = 0), as function of the 
anisotropy parameter 6. The dashed line represents the en- 
ergy of the classical model, the solid line includes corrections 
from linear spin wave theory. The extrapolated values from 
finite-size scaling of the ED results are shown as dots. Dashes 
denote the result Eq — Jc In 2 from Bethe ansatz calculations 
for the one-dimensional 5 = 1/2 spin chain. 



The dashed line shows the classical energy, the solid 
line displays the result from linear spin wave theory 
which includes corrections due to zero-point fluctuations 
of magnonsSi. The dots denote our scaled ED results 
Eo = Eo{(l},9) according to Eq. (H^ . 

Inside the magnetic phases, the overall agreement be- 
tween the numerical data and the results obtained from 
linear spin-wave calculations^'^ (solid line) is surprisingly 
good, and improves on the comparison with exact di- 
agonalization results from just a single cluster we have 
made previously^!. However, in the disordered regions 
at the borders of the columnar phase, which is shaded 
with gray, the fit to the numerical data has a compara- 
tively poor quality, and the reliability of the numerical 
result is somewhat questionable. Linear spin-wave the- 
ory breaks down here, too, albeit in a slightly different 
parameter rangeSi. 

Next we turn to the unfrustrated (J2 — 0) but 
anisotropic (Jia 7^ Jib) model with only next-nearest 
neighbor interactions. Fig. [T3] shows the dependence of 
the ground-state energy per site Eq on the anisotropy pa- 
rameter 9 for fixed frustration angle = 0. The classical 
ground-state energy is shown as a dashed line for com- 
parison. As before, the agreement with linear spin-wave 
theory (solid line) inside the magnetic phases is good, 
which is not the case at the borders of the Neel phase. 

But in this case, this has a reason different from the 
frustration induced by competing interactions discussed 
before: For 9 — or 6 — tt/2, the nearest-neighbor ex- 
change along one particular spatial direction vanishes, 
Jio or Jif,, respectively, such that we are actually dealing 
with an array of decoupled spin chains instead of a two- 
dimensional system. The classical approximation fails to 
describe the ground-state of the chains completely, and 
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FIG. 15. Structure factor Sn{Q) for tile 28:2-4 (inset: illus- 
tration of the compact tile) at the wave vectors correspond- 
ing to the four classically ordered phases as a function of the 
frustration angle (j) for the isotropic model. The labeling is 
explained in the inset of the figure. 

LSW corrections to it due to zero-point fluctuations are 
largest at these two points, improving the classical result 
drastically. The extrapolated ground-state energy from 
ED gives an even better approximation to the true value 
for Eq = Jc In 2 derived from Bethe ansatz results^, dis- 
played as short horizontal dashes in Fig. [Ml but also 
here the agreement is limited. However this is an ex- 
treme case, and we discuss it here primarily to show the 
limits of the finite-size scaling method when applied to 
the strongly anisotropic model. 

B. Structure factor and ordered moment 

The calculated values for the structure factor for all 
four ordering wave vectors as obtained using the Eq. (|24p , 
for the largest system we are considering (tile 28:2-4 with 
28 spins) is shown in Fig.tTSl Starting from the FM phase 
with a fully polarized ground state, we have Sn{Qfm) = 
S"^, and at the antiferromagnetic wave vectors Q we get 
Sn{Q) = 1/60 in accordance with Eq. 

Since the shape of tile 28:2-4 is a parallelogram, and 
hence has C2 point-group symmetry only, the equiva- 
lence between the two wave vectors Q = (tt, 0) and (0, tt) 
present for the infinite system (or any tile having at least 
C4 point-group symmetry) is lost. This manifests itself 
in a different dependence of the two structure fac- 
tors ^(QcAFa) 7^ 'S'(QcAFb) in the columnar phase, see 

Moreover, the discontinuity between the value of the 
structure factor at the (0,0) ordering in the FM phase 
and the (tt, tt) ordering in the NAF phase {(j> w — 7r/2) 
is a finite-size effect and will be suppressed by increasing 
the cluster size. 

In Fig. [16] the extrapolated values for the ordered mo- 
ment NP{Q) are plotted (with dots) as a function of the 
frustration angle <f> at two different anisotropy parame- 
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ein = 0.25 
X.. NAF 




FIG. 16. The extrapolated ordered moment as function of 
the frustration angle, for (top) the isotropic 9 = n /A case and 
(bottom) the maximally anisotropic case with 9 = Q. Those 
tiles with the maximum squareness are used for the scaling. 
The gray-shaded areas in the top plot represent the range 
of frustration angles 4> where the relative error of M^{Q) is 
above 0.1. 



ters. The top (bottom) plot corresponds to the isotropic 
(maximally anisotropic) case. Again the overall agree- 
ment with linear spin- wave theory (solid lines in Fig. [TO)) 
is very good inside the ordered regimes of the phase di- 
agram, and there are differences around the borders of 
the columnar phases. The gray-shaded areas in the up- 
per plot denote the range of frustration angles (/) with a 
relative error Croi > 0.1 (see Eq. (01])), indicating that in 
the magnetically disordered regions the numerical results 
tend to become unreliable. 

In the spin-nematic region of the isotropic model, a 
qualitative difference exists between linear spin wave the- 
ory and exact diagonalization: With increasing (/), the 
former yields a tiny region around the classical CAF/FM 
phase boundary where the ordered moment vanishes be- 
fore jumping to saturation in the FM phase. In contrast, 
the extrapolated values for M{Q) from our ED calcula- 
tions remain finite at any frustration angle 4>. 

This behavior is displayed with greater detail in 
Fig. [T71 which shows M {Q) as a function of 4> around 
the classical phase boundary, both for linear spin-wave 
theory (solid line) and exact diagonalization (dots). The 
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FIG. 17. Ordered moment M(Q) as a function of the frustra- 
tion angle, for the isotropic 9 = -r/A: case around the nematic 
region between columnar and FM phases. The solid line indi- 
cates the result from linear spin wave theory, the dots display 
the scaled ED values. The gray area shows the region where 
the relative error of the extrapolated value for M^{Q) is above 
0.1. 



extrapolated moment shows a nearly constant cj) de- 
pendence M((5cAFa) ~ M'^=°{Q^af) in the columnar 
phase for ^/tt < 0.855, as it does in the FM phase 
with M((5fm) = 5* for ^/tt > 0.874. However, for 
0.855 < (t)/7r < 0.874, Mn{Q) suddenly shows erratic 
behavior for different tile sizes (Fig. [TT] shows this for the 
long-distance correlation function), correspondingly we 
have eroi > 0.1 (no significant digits for Af^((5)) in that 
range of the frustration angle. 

The point — 0.8747r, where we can extrapolate 
Mpf{Q) to a stable limit (full polarization) again can be 
regarded as an upper bound of the border between spin 
nematic and FM phase. According to Fig. [TI] the min- 
imum value for where we get a stable M{Q) — S de- 
creases as a function of tile size, disregarding the smallest 
eight-site tile. However the true order parameter for the 
spin nematic phase is not of magnetic type, and there- 
fore the behavior of M(Q) cannot be used to understand 
the properties of this phase, in particular the parameter 
range within which it exists. Instead, one would have to 
calculate the spin nematic order parameter^- in a similar 
way, which is beyond the scope of the present paper. 

Finally, the dependence of the ordered moment on the 
anisotropy angle 9 for the unfrustrated case with only 
nearest neighbor interaction (0 = 0, equivalent to J2 = 0) 
is shown in Fig. [THl The dots display the ED extrapola- 
tion, the solid line shows the linear spin wave result. At 
the values — 7r/4 and 6 = — 3/47r, the isotropic model 
is recovered with the well-known values M((5naf) ~ 0.3 
and M((5fm) = S = 1/2. We point out again that the 
points 9 — Q and 9 — n/ 2 at the borders of the Neel phase 
correspond to arrays of independent chains. Therefore 
the moment suppression at these points is not a frustra- 
tion effect but a result of the effective one-dimensional 
character of the model, where no magnetic moment ex- 
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FIG. 18. Ordered moment M{Q) for fixed frustration an- 
gle <j) — 0, equivalent to J2 = 0, i.e., vanishing next-nearest 
neiglibor exchange. Control parameter here is the anisotropy 
angle 6; the points 9 — 7r/4 and 9 = — 3/47r denote the 
isotropic case. The solid line denotes the result from linear 
spin-wave theory, the dots display the values extrapolated 
from our ED calculations. 

ists at any wave vector. 

VI. SUMMARY 

In this work we have investigated the anisotropic frus- 
trated square lattice Heisenberg model using the exact 
diagonahzation method for finite clusters applying a fi- 
nite size scaling procedure. The latter is essential to ob- 
tain reliable values of ground state energy and ordered 
moment size in the thermodynamic limit. We have also 
compared the numerical results with those of linear spin 
wave theory^. The model which we have considered is 
characterized by an anisotropy parameter for the near- 
est neighbor interaction and a frustration angle (j) char- 
acterizing the ratio of next nearest and nearest neighbor 
exchange. 

The implementation of a stable finite size scaling pro- 
cedure requires precise criteria for the usefulness of the 
many possible clusters of varying size and shape used 
to tile the lattice. We have introduced and described 
in great detail how all possible tilings with a given area 
N of the square or rectangular lattice can be generated. 
Then we have classified the clusters according to their 
spatial symmetry, compatibility with classical magnetic 
phases and their geometrical compactness or squareness. 
We have found that a restriction to tiles which have com- 
patibility with classical phases and maximal squareness 
lead to a very stable scaling behavior of ground state en- 



ergy and ordered moment in the region of NAF and CAF 
phases. Close to the classical NAF/CAF and CAF/FM 
boundaries shown in Fig. [T] the relative error of energy 
and moment increases and the scaling procedure becomes 
unstable, because a systematic dependence of Eqn and 
Mn{Q) on the tile area N ceases to exist. In these re- 
gions, a quantitative prediction of the size of the ordered 
moment becomes difficult, if not impossible. The frus- 
tration effects of competing exchange interactions lead 
to large quantum fiuctuations which in turn cause the 
breakdown of magnetic order. 

Using the scaling results we were able to calculate the 
systematic dependence of ground state energy and or- 
dered moment as function of frustration and anisotropy 
angles (p and 9 respectively. We found that in the colum- 
nar phases, introducing a spatial anisotropy strongly sta- 
bilizes the ordered moment. In fact it becomes larger 
than that of the conventional unfrustrated isotropic near- 
est neighbor NAF with M{Q) w 0.3. This stabilization 
effect becomes particularly pronounced in the CAF/FM 
transition region where the spin nematic phase of the 
isotropic model has been found. 

The agreement of exact diagonahzation results with 
spin wave calculations was found to be generally good. 
Both methods predict the breakdown of magnetic order 
in the transition regions at the borders of the colum- 
nar magnetic phases as function of frustration but also 
in the regions where the model attains effective quasi- 
one-dimensional character as function of the anisotropy 
for J2 — 0. As in earlier investigations for the isotropic 
model, it remains difficult to quantify the exact size of 
the interval on the or 9 axes where the ordered mo- 
ment vanishes. Although not discussed here, the present 
work also has given a clear framework in which the field 
dependence of the ordered moment may be discussed. 



Appendix A: All distinct eight-site tilings of the 
square lattice 

To demonstrate how to find all possible tilings with 
a given area N of the square lattice. Fig. [1^] displays a 
list of all the distinct lattice tilings for two-dimensional 
tiles of size N = 8 created applying Eqs. ^ and dH). In 
each row, the HNF matrix H is shown together with the 
HNF tile Tff, the compact tile Tmc and the unimodu- 
lar matrix needed to map H onto M^. Note that 
in order to discuss orthorhombic or trigonal symmetry 
breaking of the Hamiltonian, we do not regard tiles as 
identical which are related by a point-group operation 
on the square lattice. 
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FIG. 19. All 15 possible distinct square lattice tilings A^f for a given tile size N = 8. The columns display the HNF matrix 
H, a graphical sketch of the HNF tile, a graphical sketch of the corresponding compact tile, and the unimodular matrix 
mapping H onto Mc. 
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